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CN Abstract: Features in the primordial power spectrum have been suggested as an 

00 explanation for glitches in the angular power spectrum of temperature anisotropies 

measured by the WMAP satellite. However, these glitches might just as well be 
CN| artifacts of noise or cosmic variance. Using the effective A% 2 between the best-fit 

CN power-law spectrum and a deconvolved primordial spectrum as a measure of "feature- 

ness" of the data, we perform a full Monte-Carlo analysis to address the question of 
how significant the recovered features are. We find that in 26% of the simulated data 
sets the reconstructed spectrum yields a greater improvement in the likelihood than 
^ for the actually observed data. While features cannot be categorically ruled out by 

this analysis, and the possibility remains that simple theoretical models which pre- 
dict some of the observed features might stand up to rigorous statistical testing, our 
results suggest that WMAP data are consistent with the assumption of a featureless 
power-law primordial spectrum. 
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1. Introduction 

The advent of a multitude of cosmological precision data in the past decade has led to 
the emergence of the so-called concordance, or "vanilla" model of cosmology. A key 
part of the vanilla model is the assumption that the spectrum of primordial curvature 
perturbations is smooth and featureless, and can be described by a simple power-law 
- consistent with their origin from an earlier period of slow- roll inflation [1-4]. 

Nonetheless, the inflationary mechanism also allows for more complex shapes 
of the spectrum, caused for instance by non-standard initial conditions [5-8], or 
unusual dynamics of the inflaton field due to, e.g., a phase transition [9,10], non- 
smoothness of the inflaton potential [11-13] or particle production [14,15] during 
inflation (see also [16-22]). In any case, any detection of features in the spectrum, 
i.e., a deviation from the standard power-law behaviour, would yield invaluable clues 
on the physics of the early Universe. Presently, the most powerful single source 
of information about the primordial state of perturbations are the observations of 
the temperature anisotropies of the Cosmic Microwave Background (CMB) by the 
WMAP satellite [23-25]. Since the first data release, a number of efforts have been 
undertaken to ascertain the compatibility of the data with the power-law paradigm. 

In a top-down approach, specific non-smooth models of the primordial spec- 
trum or the inflaton potential have been fit to the data [26-39]. Alternatively, the 
bottom-up approach of trying to reconstruct the shape of the primordial spectrum 
from these data has been employed, involving for example binning of the primordial 
spectrum [40-42], principal component analysis [43], or a direct reconstruction via 
deconvolution methods [44-53]. 

Generally, these results have shown that by introducing suitable features in the 
primordial spectrum, the fit to the data can be improved by Axis ~ 0(10) over 
the power-law fit. Taken by itself, however, A%g ff does not really help answer the 
crucial question whether the better fit indicates a real feature in the primordial 
spectrum or whether it just stems from over-fitting the scatter in the data from 
noise and cosmic variance. Indeed, attempts at interpreting results of parameter 
inference in terms of 'ruling out the power-law model in favour of a feature model' are 
likely to be prone to grossly overestimating the 'evidence' for features. A simplistic 
goodness-of-fit x 2 /dof analysis would not be very enlightening either, due to the 
large number of degrees of freedom involved. One might thus be tempted to turn 
to the formalism of Bayesian model selection [54, 55] . Though aside of the potential 
technical problems of evaluating Bayes factors for feature models with a large number 
of free parameters, the interpretation of results from a Bayesian model selection 
analysis may be problematic [56], owing to the absence of a unique well-motivated 
choice of priors for these mostly empirical parameterisations. 

We believe that for this particular problem the toolbox of frequentist statistics 
provides a useful implement, and propose to perform an analysis based on the tech- 
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nique of hypothesis testing. Given the null hypothesis of an underlying power-law 
spectrum, this involves the evaluation of a suitably chosen statistic on a large num- 
ber of Monte Carlo simulated mock CMB temperature anisotropy data sets. From 
the resulting frequency distribution of this statistic one can derive a p- value, i.e., 
the probability that, given the null hypothesis, the value of the statistic is larger 
than the one observed. We note that in defining the statistic one has to beware of a 
posteriori interpretations of the data; a particular feature observed in the real data 
may be very unlikely (and lead to a low p- value; see e.g., [57, 58]), but the probability 
of observing some feature may be quite large. 

In Section 2 (and in the Appendix) we describe in detail the methods used in our 
analysis, the results of which are presented in Section 3. We discuss and interpret 
our findings in Section 4. 

2. Method 

2.1 Null hypothesis 

We assume as our null hypothesis that the temperature angular power spectrum in- 
ferred from WMAP can be explained by the present cosmological standard ( "vanilla" ) 
model. In particular, the underlying primordial spectrum of curvature perturbations 
is taken to have a smooth power-law form 

V n (k) = A s (k/K) ns -\ (2.1) 

determined by the two parameters ng, the scalar spectral index, and As, the am- 
plitude of fluctuations at the pivot scale fc* = 0.05 Mpc -1 . The remaining four free 
parameters of the vanilla model are the baryon and cold dark matter densities 
and u c , the ratio of sound horizon to angular diameter distance at decoupling 9, and 
the reionisation optical depth r. 

The fiducial spectrum Cf d used to generate random realisations of WMAP5 data 
should be chosen as the maximum likelihood spectrum. However, since temperature 
data alone are not very sensitive to r, using temperature data alone leads to a value of 
r which is at odds with results of a combined temperature/polarisation analysis. We 
therefore fix r to the best-fit result for the full WMAP5 likelihood function, and then 
determine the values of the other five parameters by fitting the thus reduced five- 
parameter vanilla model to the temperature data only (using the likelihood function 
described in Appendix A.l). 

The resulting fiducial spectrum is defined by the following parameter values: 
co h = 0.0224, co c = 0.109, 9 = 1.04, r = 0.089, A s = 2.147 • 10^ 9 and n s = 0.963. 

2.1.1 Alternative hypothesis 

The null hypothesis is to be tested against an alternative hypothesis. Our alternative 
hypothesis is that the primordial spectrum is not given by a power-law, but possesses 
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features of some sort. It now remains to define a suitable test statistic that can be 
used to assess the validity of the null hypothesis. 

2.2 Statistic 

By relaxing the assumption of a precise power-law form of the primordial spectrum 
one can generally achieve a better fit to the data. A commonly used measure of this 
improvement is the effective delta-chi-squared 

A Xc 2 ff = -21n£ max + 21n£l x , (2.2) 

where £„L„ is the maximum likelihood of a fit to the vanilla model, and £^L Y the 
maximum likelihood of a fit to a model with an alternative shape of the spectrum. 

In the present work we wish to allow a very general form of the primordial spec- 
trum, so we focus on a reconstruction of the primordial spectrum from a binned 
version of the temperature angular power spectrum with the aid of a modified 
Richardson-Lucy (RL) deconvolution algorithm as performed in Refs. [45, 48] (see 
Appendix A. 2). Given a transfer function Te(k) and an observed or simulated angu- 
lar power spectrum Cg, this algorithm results in an "optimised" primordial spectrum 
V% L (k), and is a very powerful tool to find potential features in the spectrum. We 
take £max to be the likelihood of the data for the primordial spectrum Vn L (k) com- 
bined with the best-fit transfer function of the vanilla model. 

This method does have two limitations though: as long as the transfer func- 
tion is held fixed, one is only sensitive to features that are not degenerate with the 
cosmological parameters that determine the transfer function (e.g., typically sharp 
features such as spikes or oscillations, but not broad distortions). It would be possi- 
ble to circumvent this problem by also optimising the transfer function, but such a 
procedure would become computationally prohibitively expensive in the context of 
our analysis. Also, since we deconvolve binned angular power spectrum data, we are 
not sensitive to extremely high-frequency modulations of the primordial spectrum 
as reported in Ref. [39]. Nevertheless, the Axis obtained with our method provides 
a good measure of the potential deviation from smoothness for realistic primordial 
spectra. 

2.3 Numerical implementation 

Applied to real data, the calculation of Axlg involves a maximisation of the likelihood 
using the reduced (i.e., r = 0.089) vanilla model, and the application of the RL- 
deconvolution algorithm with the best-fit parameter transfer function. 

Naturally, the simulated data sets should be treated in exactly the same way 
as the real data. Schematically, the procedure to be performed for each random 
realisation of WMAP data should hence read as follows: 

(i) Generate a random realisation ^Cf m of simulated WMAP temperature data. 
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(ii) Find minimal xlsv f° r a fit of the vanilla model to wC| im . 

(iii) Calculate the transfer function ®Ti(k) for the best-fit parameters found in step 
(ii). 

(iv) Apply the modified Richardson-Lucy deconvolution algorithm to ^'Cf m , using 
the transfer function ( l 'Tg{k). 

(v) Determine xls rl fr° m the deconvolved primordial spectrum and calculate the 
improvement A%g ff over the power-law fit. 

To keep the task of simulating mock WMAP data at the level of generating angular 
power spectra ^Cf™ and avoid having to generate random maps, we adopt a slightly 
simplified form of the WMAP likelihood function (described in Appendix A.l) through- 
out this paper. For details of simulating mock WMAP data (step (i)) we refer the 
reader to Appendix A. 3. 

The minimisation of step (ii) turns out to be the computationally most expensive 
part of the algorithm. Deterministic, simplex-like minimisation routines like Numer- 
ical Recipes' amoeba [59] are not very reliable for our purpose since they run the 
risk of getting stuck in local minima. Random-walk based routines, such as simu- 
lated annealing are more suitable here. We choose a combination of the amebsa [59] 
routine to find a good starting point for a subsequent low-temperature run of the 
Metropolis-Hastings algorithm (based on the CosmoMC code [60]), which results in 
an estimate of the minimal Xcfrv that is accurate to about 0.1. Nonetheless, this 
algorithm requires several thousand evaluations of CMB angular power spectra and 
the corresponding likelihoods, which would take a considerable time if one were to 
use a conventional Boltzmann code, such as CAMB [61]. We resort instead to the 
interpolation code PICO [62, 63] for calculating the Ces, which can considerably speed 
up this task, returning a reliable estimate for the minimal xis v within a few minutes 
on a regular desktop CPU. 

Having found a good estimate of the maximum likelihood point, we then use 
CAMB to calculate the transfer function for the best-fit cosmological parameters in 
step (iii). Having applied the RL-deconvolution (step (iv)), giving us the optimised 
spectrum V-r, (k), we can finally compute the improvement ^xls the optimised 
spectrum would yield over the best power-law one (step (v)). 

3. Results 
3.1 Real data 

We commence by applying the RL-deconvolution method to the observed WMAP5 
temperature data. The resulting best fit "optimised" primordial spectrum is shown 
in Fig. 1; its dominant features are qualitatively similar to those found in Ref. [48] for 
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the 3-year WMAP data: a cutoff-like suppression at the largest scales accompanied by 
two prominent wiggles around k ~ 0.002 Mpc -1 and k ~ 0.07 Mpc -1 . The optimised 
spectrum gives xIsrl = 1029.45, which marks an improvement of A% 2 = 24.41 over 
the power-law best fit of xts v = 1053.86. 
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Figure 1: Best-fit power-law (black line), and RL-optimised primordial spectrum (red 
line) . 

3.2 Simulated data 

We have generated 2000 random realisations of simulated WMAP5 temperature data 
for multipoles 2 < £ < 1000 and applied the procedure of Section 2.3 to determine 
the A%g ff for each of them. A histogram of the distribution of A%g ff values is shown 
in Figure 2. Since the RL-algorithm is applied to binned data whereas Axg ff is 
defined using un-binned data, for some simulated spectra the resulting A%g ff can 
be negative, even after wavelet smoothing. In other words, in these cases the RL- 
optimised spectra yield no improvement over a power-law, and we set A^g = 0. 
This has occurred in about 5% of our simulated spectra and is marked by the red 
part of the lowest-Ax^-bin in the histogram. 

The vertical blue line in Figure 2 marks the observed value Ax^s = 24.41. Of the 
2000 mock data sets, 525 have A%g ff > 24.41, corresponding to a p-value of ~ 0.26. 
To ensure this number is not affected by our simplified likelihood-approximation, we 
calculated the observed Ax^f using the full WMAP5 likelihood code (version v3p2), 
and found a difference of ~ 1. 

Additionally, we considered statistics which are less susceptible to low multipoles 
where the likelihood approximation is most inaccurate: Ax^ ff evaluated between £ min 
and £ max = 1000 with £ min = 10, 20, and 30. The resulting p-values are 0.33, 0.14, 
and 0.31, respectively. This behaviour can be understood from the features observed 
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Figure 2: Histogram showing the relative frequency of ^-Xcs f° r 2000 simulated realisa- 
tions of the WMAP5 temperature spectrum. For about 5% of the spectra, the RL-algorithm 
does not result in a better fit; these cases are displayed in red. The dotted blue line denotes 
the observed value Ax^g- = 24.41. 

in the real data: removing i < 10 gets rid of a feature, leading to an increase in 
the p-value; taking away 10 < i < 20 removes a part of the data which does not 
show features, causing the p-value to drop, only to go up again when one removes 
the feature at 20 < I < 30. 

In summary, the improvement in the fit from using an "optimised" primordial 
spectrum over a power-law spectrum is larger than the observed one in roughly 30% 
of our simulated data sets which assume an underlying smooth spectrum. We thus 
conclude that the null hypothesis cannot be rejected at high significance within the 
limitations of our chosen statistic, and the data are compatible with a smooth power- 
law primordial spectrum. We illustrate this qualitatively in Fig. 3, which shows that 
the reconstructed spectrum from the observed data does not have an unusual amount 
of "featureness" compared to the simulated data. 

4. Conclusions 

We have addressed the question whether the 5-year WMAP temperature anisotropy 
data are compatible with the assumption of an underlying smooth power-law pri- 
mordial spectrum of curvature perturbations, or whether they show any indication 
for features or other unaccounted-for systematic effects. Assuming the true under- 
lying primordial spectrum to be of power-law form, we generated 2000 simulated 
WMAP angular temperature power spectra, and estimated the amount of features 
by looking at the improvement to the likelihood gained from fitting an "optimised" 
primordial spectrum, obtained by deconvolving the data, instead of the usual power- 
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Figure 3: Reconstructed primordial spectra for twenty mock data sets (orange lines). 
For comparison we show the primordial spectrum reconstructed from WMAP5 data (red 
line) and our fiducial power-law spectrum (black line). It is evident that the spectrum 
reconstructed from real data does not have an unusual amount of features. The apparent 
feature at 0.05 Mpc -1 < k < 0.07 Mpc -1 is caused by the noise term becoming dominant at 
the corresponding multipoles in the binned WMAP data. For the last bin, 951 < I < 1000, 
the averaged transfer function peaks at k ~ 0.07 Mpc -1 . Beyond that the reconstructed 
primordial spectrum is no longer dominated by the peak of one data bin, but rather by 
the tails of the transfer functions of the last few bins; more structure here does not further 
improve the fit. 

law. We found that 26% of all simulated data sets show a greater improvement than 
the one observed, which leads us to conclude that the features seen in the WMAP 
5-year temperature data are not at all incompatible with the assumption of a smooth 
primordial spectrum. 1 

However, we emphasise that our analysis does not disprove the existence of 
features - they are merely not strictly required by present data. Specifically, 'simple' 
theoretical models that predict similar features could conceivably remain favoured 
over a power-law. Additionally, since the spectrum reconstruction method used here 
is not sensitive to features below the binning scale of the WMAP data, we cannot 
rule out the possibility of extremely high-frequency features. We have demonstrated 
the feasibility of a full-scale frequentist hypothesis testing analysis in the search for 
deviations from smooth spectra though, and our method can just as well be applied to 
other statistics which allow for even more general shapes of the primordial spectrum. 

Future analyses and new measurements may well reveal evidence for the existence 

1 While this work was being completed, Ref. [65] appeared, which reaches similar conclusions 
from cross- validating a spline-reconstruction of the primordial spectrum from various cosmological 
data sets. 



-7 - 



of features, and a search for them will certainly remain a worthwhile endeavour, 
considering that any detection would profoundly impact our understanding of the 
physics of inflation and also precise estimation of the cosmological parameters [64]. 

Besides the imminent improvements to temperature anisotropy data from the 
Planck satellite [66], better measurements of the .^-polarisation of the CMB will 
greatly enhance sensitivity to primordial features [52, 67] - a dedicated mission like 
CMBPol [68] may prove particularly helpful here. 
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A. Appendix: Technical details 
A.l Likelihood function 

For reasons of simplicity we adopt here a form of the likelihood function similar to the 
one used by the WMAP team in their first data release [69] . For an input theoretical 
angular power spectrum Cf 1 , the likelihood of the data reads 




(A.l) 



with a Gaussian part 




(A.2) 
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and a log-normal part 



-2\nC' LN = J2(4 h - Ze){Cf+Me)Qee>{Cf+M v ){zf - z e ), (A.3) 
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where zf 1 = \n(Cl h + Me) and ze = ki(Ce + Me), with Ce the angular power spectrum 
estimated from observation. The curvature matrix Qw is given by 



Qee = L>e <w + 7= = , 
V LftlJt' 
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with 

°' ~ 2 WTT)Ji r - (A ' 5) 

The off-diagonal terms induced by the sky cut r«/, the effective noise power spectrum 
Afg and the effective sky fraction / S k y are supplied by the WMAP team [23] and 
available for download on the lambda web site 2 . 

A. 2 RL-algorithm for reconstructing the primordial spectrum 

The Richardson-Lucy (RL) algorithm was developed and is widely used in the context 
of image reconstruction in astronomy [70,71]. However, the method has also been 
successfully used in cosmology, to deproject the 3-D correlation function and power 
spectrum from the measured 2-D angular correlation and 2-D power spectrum [72, 
73]. 

The angular power spectrum, Cg, is a convolution of the initial power spectrum 
V-ji(k) generated in the early universe, with a radiative transfer kernel Tg(k) (the 
transfer function), that is determined by the values of the cosmological parameters. 
In our application, we solve the inverse problem of determining the primordial power 
spectrum, V-ji(k), from the measured angular power spectrum, Cg, using the relation 

Ci = f dkT £ (k)V n (k) ~ ^T^P^ki). (A.6) 

In the above equation, the target measured angular power spectrum, Cg = Cf° s , 
is the data given by observations, and the radiative transport kernel, 

Ah- 

Tt{h) = — ^|A T ,(/c i3 r7o)| 2 , (A.7) 

encodes the response of the present multipoles of the CMB perturbed photon distri- 
bution function A^/cj,^) to unit of power per logarithm interval of wavenumber, 
k, in the primordial perturbation spectrum. The kernel Tg(k) is completely fixed by 
the cosmological parameters of the base cosmological model. The kernel Tg(k) also 
includes the effect of geometrical projection from the three-dimensional wavenumber, 
k, to the harmonic multipole, I on the two-dimensional sphere. 

Obtaining V-jz(k) from the measured Cg for a given Tg(k) is clearly a deconvolution 
problem. An important feature of the problem is that C| , Tg(k) and Vn(k) are all 
positive definite. We employ an improved RL method to solve the inverse problem 
for Vfz{k) in Eq. (A.6). The advantage of the RL method is that positivity of the 
recovered P^(fc) is automatically ensured, given positive definite Tg{k) and Cgs. The 
RL method, readily derived from elementary probability theory of distributions [70] , 
is an iterative method that can be neatly encoded into a simple recurrence relation. 
The power spectrum V-r,^ (k) recovered after iteration (i + 1) is given by 
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where Tg{k) is the normalised kernel (on the i space for all k wavenumbers), C° bs 



(Ai 

^o: 

'I 



is the normalised measured data (target) and Cjfi is the angular power spectrum at 
the i th iteration obtained from Cf = J2T e (k)V n {i) (k) using the recovered power 
spectrum V-n (k). It is important to keep in mind that, due to the formulation 
in terms of conditional probability distributions, the RL method requires the kernel 
Ti(k), data Ce, and the target vector Vn{k), all to be normalised at the beginning, 

53 C e = 1; V ~k {1) W = E f ^ = ( A ' 9 ) 

where Vn l \k) is the normalised initial guess model of the primordial spectrum and 
the normalisation factor is defined by = Yle Te{k) . 

Note that Eq.(A.8) represents a modified form of the Richardson-Lucy algorithm 
specifically tailored for our purpose in this problem [45,47,48]. 

Following Refs. [45, 47, 48], we apply the RL-algorithm to a binned version of the 
(real and simulated) WMAP5 data. The resulting raw deconvolved spectrum has 
spurious oscillations and features on scales smaller than the bin size, arising largely 
due to the k space sampling and binning in i space, which adversely impact the 
(unbinned) full likelihood of the data given the reconstructed spectrum. As shown in 
Ref . [47] , a subsequent smoothing of the raw reconstructed spectrum with a discrete 
wavelet transform can lead to a significant improvement in the likelihood. 

A. 2.1 Discrete Wavelet Transform 

Wavelet transforms provide a powerful tool for the analysis of transient and non- 
stationary data and is particularly useful in picking out characteristic variations at 
different resolutions or scales. This linear transform separates a data set in the form of 
low-pass or average coefficients, which reflect the average behaviour of the data, and 
wavelet or high-pass coefficients at different levels, which capture the variations at 
corresponding scales. As compared to Fourier or window Fourier transform, wavelets 
allow optimal "time-frequency" localisation in the real as well as in the Fourier 
domain. The vocabulary of DWT stems from applications in one dimensional time- 
stream signal trains, but has found wide application in signal in other domains and 
dimensions. Specifically in our case, the signal being transformed is the primordial 
power spectrum, Vji(k), a one dimensional function of wavenumber, k. 

Wavelets are an orthonormal basis of small waves, with their variations primar- 
ily concentrated in a finite region, which make them ideal for analysing localised 
transient signals. Wavelets can be continuous or discrete. In the latter case, the 
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basis elements are strictly finite in size, enabling them to achieve localisation, while 
disentangling characteristic variations at different frequencies [74]. For more details 
about DWT and its theoretical basis, see [47]. 

In this paper, we use the prescription of Ref. [48], and smooth the raw decon- 
volved spectrum with the following method: 

1. Increase the number of A;- values of the spectrum to a power of 2, by padding 
evenly on both sides 

2. Perform discrete wavelet transform 

3. Keep the 2 N wavelet coefficients corresponding to lowest frequency features 

4. Perform inverse discrete wavelet transform and calculate the likelihood for the 
smoothed spectrum 

5. Repeat steps 3.-4. for N G {0,11} and keep the spectrum that provides the 
best likelihood. 

A. 3 Generating random realisations of WMAP data 

Given an observed CMB temperature anisotropy AT(n)/T, after performing the 
usual expansion into spherical harmonics 

AT(h)/f = J2&m(n), (A.10) 

£m 

we can define an observed angular power spectrum 

1 1 

m=— I 

which can be split up into an original signal C° bs and a contribution from experimental 
noise A/^ obs . Theory, on the other hand, can only predict the average of this quantity 
over an ensemble of independent observations, Cf 1 = (<£'f' s )- 

In the following, we shall describe how to generate independent realisations of 
simulated mock data (which replace the Ce in the associated mock likelihood function) 
from a theoretical input spectrum (£f d . Since the Ce are actually defined as the 
difference between the observed spectrum and the expected noise power spectrum, 
Ce = ^£ bs — A/^ th , the corresponding simulated quantity is <tf m — Mf 1 . 

For the sake of clarity we commence with the case of an ideal (i.e., noise-free, 
full-sky) observation. 



- 11 - 



A. 3.1 No experimental noise, full sky-coverage, no correlations 

Even in the case of an ideal measurement of the CMB temperature anisotropies 
(Me = 0), the observed angular power spectrum would still be subject to cosmic 
variance. Under the assumption of Gaussian fluctuations, the coefficients d£ m can be 
considered Gaussian random variables for the purpose of simulating data sets Cf m . 
Isotropy dictates (a£ m a}, m ,) = Cf d 5u>5 mm ', so different multipoles will be independent 
of each other. Hence, for each i the simulated C| im can be constructed by drawing 
random numbers xg from a xf^+i-distribution: 

A. 3. 2 Adding experimental noise 

Realistic observations are of course subject to instrumental noise. 

Assuming the noise to be isotropic and Gaussian, its effect on the distribution 
the ai m are drawn from can be considered as a convolution with a Gaussian of width 
i/A^ th , and we have 




(A.13) 



A. 3. 3 Partial sky coverage and correlations 

Even for nominally full-sky observations, some parts of the sky cannot be used to 
construct spectra and need to be cut out, since the signal is obscured by point 
sources and galactic foregrounds. This has two consequences on the recovered C^. 
First, using only part of the sky leads to a loss of information, resulting in a larger 
scatter of Ce around the true underlying spectrum Cf d . Second, since the sky cut 
breaks isotropy, the (pseudo)-C^ constructed from an incomplete map will no longer 
be uncorrelated [75]. 

Since the curvature matrix (A. 4) appearing in the likelihood function scales like 
/ s 2 ky , the standard deviation of C| im — Cf d should scale like . This increase in 
scatter around the fiducial model can be modelled by drawing the random variables 
xt from a Xf2 ^^-distribution. Since it is numerically much easier to generate \ 2 
variates with an integer number of degrees of freedom we round / s 2 k (2f + 1) to the 
nearest integer and eventually rescale the resulting y e = xe/ (xe) — 1 accordingly, to 
recover the correct standard deviation. 

The so-generated ye are uncorrelated random variables. For the simulated spec- 
trum to have the same correlation properties as the real WMAP data, we take instead 
the vector 

{y t )=Cho\{Kur){yv\ (A. 14) 
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where the lower triangular matrix Chol(iQ^) is the Cholesky decomposition of the 
correlation matrix Kg£ = 5u> + r«/. Finally, the simulated spectra are given by 

q im - N? = C\ A + (Cf d + A/f ) y t . (A. 15) 

The above procedure can rarely result in negative values of Cf m +Nf m at low multi- 
poles, due to the rescaling of the ye or due to the transformation of Eq. (A. 15). Such 
data is unphysical and not compatible with the log-normal piece of the likelihood 
function and must therefore be discarded. As a result, the simulated data can be 
slightly biased; however, the bias is negligible and remains smaller than (9(1%) of 
cosmic variance for all £. 
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